Repair the dual basis on QUAD8 and TRI6 faces - #4536
Conversation
The weight \int N_k used to build the biorthogonal dual basis is not positive on the serendipity quadratic faces. It is exactly 0 at a TRI6 vertex, which leaves those vertex dual shape functions identically zero, and -1/3 at a QUAD8 corner, which flips the sign of any per-node normalization built from it. Neither case is currently diagnosed. Biorthogonalize instead against the locally quadratic transformed basis of Popp et al., SIAM J. Sci. Comput. 34(4):B421-B446, 2012, Sec. 4.4.1, on exactly those two element types. Each vertex absorbs a fraction alpha of its adjacent mid-edge shapes and each mid-edge shape is scaled by 1 - 2*alpha; making the dual biorthogonal to that basis replaces D by T^-1 diag(T d), which for this T reduces to a sparse update of D needing no extra linear solve. With alpha = 1/5 the transformed weights are strictly positive (QUAD8 -> 1/5 and 4/5, TRI6 -> 1/15 and 1/10) and the partition of unity is preserved, so the dual basis still reproduces constants. The transform is the identity on every other element type, and the alpha = 0 limit reproduces diag(d) exactly, so all other dual bases are unchanged. dual_shape_verification_test gains QUAD8 and TRI6 cases covering positivity, reproduction of constants, and biorthogonality against the transformed basis. Refs idaholab/moose#33558 Co-authored-by: Claude <noreply@anthropic.com>
|
I've started up CI. Walk me through the math on this, and make sure I'm not misunderstanding? The standard dual basis on these elements can't be made to satisfy both biorthogonality and the partition-of-unity property at once without going degenerate, so we do a slight change-of-basis on the primal bases and then we find a dual to that changed basis instead? How does that work at the application level, when we no longer have a dual to the original basis but we aren't exposing the altered basis? I think I'm okay with this, under the circumstances (though it pains me to see a non-serendipity element behaving as badly and requiring as much of a hack as a serendipity!), but we've got to document what's going on in the headers, not just the implementation. Somehow I was an idiot and approved the original addition of |
|
Partition of unity isn't the obstruction. In §4.4.1 of (Popp, Wohlmuth, Gee, Wall) the authors note the transformation is chosen to preserve it. What fails is their integral positivity condition, Eq. (4.2): ∫Φ_j dγ > 0. Eq. (4.3) gives < 0 at a quad8 corner and Eq. (4.4) gives exactly 0 at a tri6 corner, where they say biorthogonality "does not even make sense". So tri6 is the genuinely degenerate case; quad8 just has the wrong sign. |
get_dual_phi() and its siblings carried no Doxygen at all, and on QUAD8 and TRI6 they now return a basis biorthogonal to a transformed primal basis rather than to N, so \int Phi_j N_k is not diagonal on those two types. Spell that out where users will look, including that the mortar matrix still factors as D = D~ T^-1 with D~ diagonal, leaving D^-1 = T D~^-1 sparse and cheap to apply. Also stop calling TRI6 serendipity -- it is the complete quadratic triangle -- and use the integral positivity terminology and equation numbers of the cited paper. Co-authored-by: Claude <noreply@anthropic.com>
|
Done in bc22f26.
Also dropped "serendipity" where I had wrongly applied it to tri6, and switched the implementation Full unit suite still passes locally. (generated by Claude) |
3 participants
compute_dual_shape_coeffsbuilds the dual basis fromD = diag(\int N_k), but that weight isexactly 0 at a TRI6 vertex -- leaving those vertex dual shape functions identically zero -- and -1/3
at a QUAD8 corner. Neither is diagnosed.
This biorthogonalizes against the locally quadratic transformed basis of Popp et al. (2012, Sec.
4.4.1) on those two element types, which reduces to a sparse update of
Dwith no extra solve. Thetransform is the identity everywhere else and
alpha = 0reproducesdiag(d)exactly, so nothingelse changes, there is nothing to regold, and no API changes.
New QUAD8/TRI6 cases in
dual_shape_verification_testfail without the fix. Full unit suite passes;MOOSE first-order dual mortar output is bit-identical before and after.
Refs idaholab/moose#33558